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The  report  describes,  in  some  detail,  a  proposed 
computer  programme  for  the  dynamic  analysis  of  elasto-plastic , 
thin,  isotropic  plates  of  rectangular  planfornn  Specific  dynamic 
loading  is  considered  which  is  appropriate  to  a  blastwave,  namely, 
a  spatially  constant  pressure  with  infinitesimal  rise  time, 
decaying  rapidly  to  a  suction  peak  and  finally  to  ambient 
conditions.  Account  is  taken  of  the  large  deformations  of  the 
plate,  and  fixed  or  simply-supported  boundary  conditions.  (U) 


Much  of  the  programme  has  been  coded  and  step-by-step 
testing  of  the  various  subroutines  is  currently  progressing.  A 
flow  chart  of  the  main  programme  is  given ,  iiriftppendix  I,  and  a 
statement  of  the  subroutine  testing  completed  as  at  the 
publication  date  is  given  in  Appendix  II. 
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matrix  relating  stress  coefficients  and  local  nodal 
displacements 
elasticity  matrix 
flow  matrix  defined  by 

instantaneous  slope  of  unaxial  stress  versus  plastic 
strain  plot  at  beginning  of  time  increment 
time  increment 
hybrid  stiffness  matrix 

axes  transformation  matrix  relating  to  a  node 
consistent  mass  matrix 
compliance  matrix 

total  number  of  elements  in  idealization 
number  of  Gaussian  quadrature  points  per  element 
total  number  of  nodes  in  idealization 
number  of  nodes  per  element 
vector  of  applied  nodal  loads 

matrix  relating  elastic  stress  components  and  coefficients 
vector  of  nodal  displacements,  velocities  and  accelerations, 
respectively 

vectors  used  in  accumulating  incremental  nodal  displacements 
and  velocities,  respectively 

vector  of  nodal  pseudo-forces  accounting  for  plasticity 
vector  of  accumulative  nodal  pseudo-forces  incorporating 
effects  due  to  straining  and  large  rotation  from  initial 
configuration 

axes  transformation  matrix  relating  to  an  element 
time  co-ordinate 

matrix  of  Gauss  weights  and  co-ords  over  area  of  element 
"  ”  "  "  "  "  through  thickness 

rectangular  space  co-ordin±es 
vector  of  stress  coefficients 
work  hardening  parameter 
proportionality  constant  in  flow  rule 
vector  of  stress  components 
uniaxial  stress  at  first  yield 
equivalent  stress  at  a  station 
uniaxial  yield  stress  subject  to  work  hardening 
vector  of  strain  components 
equivalent  plastic  strain  at  a  station 
uniaxial  plastic  strain 
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Subscripts 


€  elastic 

«ff  effective  (i.e.  uniaxial) 

tp  elasto-plastic 

L  relating  to  local  axes 

o  original 

p  plastic 

R  revised 

6  relating  to  additional  straining  effects 

■&  relating  to  large  rotation  effects 

Superscripts 
C  elemental 

g  at  a  Gauss  Point 

j  nodal 

n  time  step  counter 

T  matrix  transposition 
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1.  INTRODUCTION 


The  capability  of  large  modern  computers  makes  it 

possible  to  solve  physical  problems  that  defy  a  rigourous  .  .• . 

mathematical  approach.  One  such  problem  of  interest  in 
engineering  is  the  behaviour  of  thin  plates  subjected  to  blast 
loading.  Obviously,  the  analysts  of  military  systems  have 
a  keen  interest,  as  do  the  engineers  concerned  with  explosive 
forming. 

The  Finite  Element  method  of  stress  analysis  is  ideally 
suited  for  large-scale  computer  systems  and  many  successful 
analyses  of  complex  structures  have  been  reported.  Steady 
progress  from  static  to  dynamic  situations  is  apparent,  and, 
more  recently,  the  non-linear  behaviour  of  both  material 
and  structure  has  commanded  considerable  attention. 

Since  elasto-plastic  analysis  is  heavily  dependant  on 
stress  prediction,  it  is  essential  to  use  a  type  of  finite 
element  formulation  which  has  a  good  reputation  for  quite  accurate 
stress  analysis,  without  incurring  excessive  computer  time  by 
dealing  with  a  fine  element  division.  Such  a  formualtion  is 
the  "Hybrid"  technique,  developed  by  Pian>  which  falls  between 
the  strictly  displacement  methods  and  the  "stress-mode"  methods. 
Thus,  the  present  research  is  centred  on  a  well  tested 
"hybrid"  element,  which  may  be  triangular  or  quadrilateral 
in  planform. 

However,  very  considerable  adaption  of  the  element  has 
proved  necessary  due  to  revised  integration  schemes  through 
the  volume  in  order  to  account  for  the  partial  plasticity 
within  each  element.  The  5a.  Reformations  are  taken  into 
account  by  psuedo-forces  acting  at  the  nodes,  and  the  non-linear 
equations  of  motion  integrated  by  a  Runge-Kutta  predictor  - 
corrector  scheme. 


2  A  Survey  of  Solution  Methods  for  Elasto-Plasticitv  Problems 

This  survey  looks  at  recent  work  performed  on  the  solution 
of  structural  problems  involving  non-linear  material  behaviour, 
and  is  restricted  to  elasto-plasticity.  Other  material  non- 
linearities,  such  as  creep,  are  not  included  and  the  reader 
;'.s  referred  to  the  text  Ly  Zienkiewicz  (1)*  for  discussion 
of  these  phenomena. 


Firstly  we  shall  briefly  consider  some  of  the  latest  1 

experimental  work  performed  in  this  area.  A  number  of  investigators 
have  turner',  their  attention  to  beams  under  blast  loading.  For 
example,  Humphreys  (2)  has  tested  straight  clamped  beams,  and 
Florence  and  Firth  (3)  have  subjected  both  pinned  and 
clamped  rigid-plastic  beams  to  uniformly  distributed  impulses. 

Grid  frameworks  have  been  investigated  (4)  as  have  perforated 
plates  in  plane-stress  (5).  Recently  Jones  et  al.  (6,7) 


performed  tests  on  the  dynamic  behaviour  of  fully  clamped 
rectangular  plates  and  compared  their  results  (6)  with  the 
theories  of  Martin  (8)  and  Haythornthwaite  and  Shield  (9). 


Although  progress  was  being  made  in  the  mathematical 
theory  of  plasticity  (10,11)  it  was  not  until  the  advent  of 
the  high-speed  digital  computer,  and  in  particular  the 
application  of  the  Finite-Element  method  to  elasto-plastic 
materials,  that  the  analysis  of  complex  problems  of  this  type 
became  anything  but  very  approximate.  We  shall,  however, 
postpone  discussion  of  this  approach  until  we  have  considered 
some  non-Finite-Element  solutions.  Usually  major  simplifying 
assumptions  were  made  -  often  elastic  effects  were  neglected 
(12,13)  and  perfect  plasticity  assumed  (14.15).  Prager  (16) 
reviews  the  state  of  the  art  in  the  mid  1950's.  At  about 
this  time  Hopkins  and  Wang  (17)  had  analysed  circular  perfectly 
plastic  plates  and  compared  results  obtained  using  Tresca,  von 
Mises  and  Parabolic  (von  Mises)  yield  conditions.  Good 
agreement,  especially  for  the  simply-supported  condition, 
between  Tresca  and  von  Mises  was  obtained.  Ang  and  Lopez  (18) 

*  such  numbers  are  references 
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employed  a  grid  analogy  method  and  the  usual  sandwich 
plate  assumptions  in  deriving  their  plate  analysis.  Deformable 
nodes  represented  average  flexural  and  axial  resistance, 
torsional  elements  modelled  the  in-plane  sheer  stresses  and 
rigid  bars  transmitted  the  transverse  shear  forces  induced  by 
the  flexural  stresses.  Results  were  given  for  the  progression 
of  the  elasto-plastic  boundary  in  square  plates  under  various 
edge  conditions. 

Dynamic  analyses  of  elasto-plastic  plates  have  also 
been  performed  (12,  13).  In  1959  Cox  and  Morland  (19) 
considered  a  simply-supported  square  plate  subjected  to  a 
pressure  pulse,  and  later  Florence  (20)  analysed  clamped, 
circular  plates  under  a  rectangular  blast  pulse  by  assuming 
a  rigid-plastic  material  yielding  at  the  Tresca  hexagon. 
Recently  problems  of  combined  material  and  geometric  non- 
linearities  have  been  tackled.  Gerdeen  et  al. (21)  have 
analysed  shells  of  revolution  under  these  conditions  and 
Symonds  and  Jones  (22)  have  investigated  the  behaviour  of 
fully  clamped  beams  composed  of  materials  (particularly 
steels)  which  were  strain-rate  sensitive. 

The  Finite  Difference  approach  has  been  used  for  some 
time  to  analyse  elasto-plastic  problems  (e.g.  23,  24)  sometimes 
in  combination  with  large  deflections.  Balmer  and  Witmer 
(25)  performed  such  an  analysis  on  impulsively  loaded  beams, 
and,  with  high-energy  metal  forming  in  mind,  the  extension 
to  thin  shells  has  been  made  (26,  27).  Leech  et  al,  (26) 
presented  results  for  a  cylindrical  panel  and  compared  them 
with  experiments,  and  Lindberg  and  Boyd  (27)  considered 
clamped,  rigid-strain  hardening  shell  membranes.  It  is 
however,  evident  that  further  work  needs  to  be  done  in  this 
particular  field. 

So  far  the  Finite-Element  method  has  been  mentioned 
only  fleetingly.  Now  we  come  to  consider  its  application 
to  the  analysis  of  elasto-plastic  structures.  Although  the 
Force  approach  (28,  29)  has  been  used  in  this  context,  notably 
by  Denke  (30)  and  Lansing  (31),  the  vast  majority  of  workers 
have  used  the  stiffness  (displacement)  method  (32).  It  is 
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applications  of  this  approach  which  will  now  be  considered. 


One  of  the  most  basic  ways  of  solving  the  non-linear  j  \ 

equations  is  by  the  Direct  Iteration  procedure.  In  this  j  - 

method  the  total  load  is  applied  to  the  structure  and  using  I 

the  initial  values  of  Young's  Modulus,  E,  and  Poisson's  | 

RatioyVjthe  elastic  stresses  and  strains  are  computed.  \ 

\ 

New  values  of  E  and  V  are  calculated  according  to  the  stress  j 

levels  reached  in  each  element,  and  then  another  full  "elastic"  j  L 

analysis  performed  using  these  values.  This  procedure  is  j  . 

repeated  until  the  solution  has  converged.  Figure  la)  shows  { 

a  diagrammatic  representation  of  this  process  as  applied  to 

a  material  with  an  elasto-plastic  stress-strain  curve  which 

may  be  represented  adequately  by  two  straight  lines.  In 

1963  the  method  was  applied  by  Wilson  (33)  to  2-D  structures 

and  later  used  by  Clough  (34).  These  workers  demonstrated 

that  the  Direct  Iterative  approach  usually  converges 

satisfactorily  after  3  or  4  iterations.  A  completely  different 

direct  solution  has  been  demonstrated  by  Mallet  and  Schmit 

(35).  They  use  an  Energy  Search  technique  for  a  large  deflection 

analysis  and  emphasise  that  the  process  is  equally  applicable 

to  elasto-plasticity,  although  this  point  is  disputed  by 

Hofmeister  et  al.(50). 


Probably  the  first  method  used  for  Finite-Element 
elasto-plastic  analysis  was  the  incremental  Initial  Strain 
approach  introdudec^y  Mendelson  and  Mason  (36)  in  1959. 

Padlog  et  al»(37)  considered  inelastic  structures  under 
cyclic,  thermal  and  mechanical  loading  in  1960,  and  other 
pioneer  work  was  performed  by  Gallagher  et  al. (38),  Argyris 
(39)  and  Jensen  et  al. (40).  The  procedure  is  to  apply  the 
load  incrementally  and  at  each  load  step  to  treat  the  plastic 
strains  produced  as  initital  strains  for  the  next  increment. 

In  this  way  pseudo-loads  account  for  plasticity  and  it  is 


therefore  necessary  to  modify  only  the  right-hand-sides  of 

equilibrium  equations.  This  means  that  the  elastic  stiffness 
matrix  is  used  throughout  and  thus  needs  to  be  built-up  and 

inverted  once  only.  Using  the  usual  Finite-Element  notation 

we  may  write  the  linear  elastic  stress-strain  relationship 

for  a  material  as 


cr  -  crc  =  [D]  (jE  "  £?) 


In  this  approach  it  is  the  initial  strain  vector,  £3  1 
which  is  adjusted  to  achieve  a  solution  of  the  equilibrium 
equations 


=  P 


where  represents  the  vector  6f  all  forces  due  to 
external  loads,  initital  strains,  initial  : tresses  etc. 

The  procedure  is  shown  diagrammarically  in  Figure  lb).  If 
large  increments  of  load  are  taken  it  can  become  necessary 
to  iterate  at  each  step  in  order  to  achieve  high  accuracy,  but 
this  can  easily  be  avoided  by  taking  the  increments  sufficiently 
small  (39).  A  major  disadvantage  of  the  Initial  Strain  approach 
is  that  it  breaks  down  for  perfectly  plastic  materials  because 
the  strains  in  this  case  cannot  be  uniquely  determined 
for  prescribed  stress  levels. 

A  central  problem  in  an  Initial  Strain  analysis  is  the 
method  used  to  compute  the  pseudo-forces.  This  is  a  much 
more  crucial  point  in  plate— bending  than  in  2-D  analysis. 

Many  non-Finite-Element  plate-bending  analyses  assumed  that 
at  all  points  on  the  plate  the  entire  thickness  was  either 
fully  elastic  or  fully  plastic  (23,  18),  and  whilst  this  is 
the  case  for  a  plate  of  sandwich  construction  it  is  a  very 
crude  approximation  for  a  solid  plate.  Using  the  Finir'''— 
Element  method  less  drastic  assumption  can  be  made  to 
overcome  this  problem.  For  example,  Armen 'et  al.(41)  for 
their  initial  Strain  analysis  employed  a  triangular  plate¬ 
bending  element  with  quintic  displacement  functions  and 
assumed  the  plastic  strain  to  vary  linearly  from  the  upper 
and  lower  surfaces  of  the  element  to  elastic-plastic 
boundaries  within  the  cross-section.  They  further  ass^ied 
a  linear  variation  along  the  edges  between  adjacent  nodes 
with  a  corresponding  planar  distribution  throughout  the  area 


of  the  element.  Nodal  strains  were  found  by  averaging  values 
from  surrounding  elements.  In  this  work  the  Romburg-Osgood 
stress-strain .relationship  (42)  is  employed,  and  so  that  the 
Bauschinger  effect  could  be  included  for  cyclic  loading,  the 
authors  chose  the  Prager-Ziegler  kinematic  hardening  theory  of 
plasticity  (16,  43,  44).  However,  recently  this  group  of 
researchers  have  favoured  numerical  integration  techniques 
both  along  the  length  and  through  the  thickness  of  their  shell 
of  revolution  (45).  They  have  found  3  Gauss  points  generally^ 
adequate  along  the  length.  Through  the  thickness  Simpson's 
integration  employing  up  to  21*..station  was  selected  -  the 
large  number  being  required  for  the  cyclic  loading  being 
considered  in  this  work.  The  choice  of  Simpson's  integration 
allowed  the  first,  surface  yielding  to  be  detected,  whereas 
Gaussian  quadrature  would  have  required  yielding  at  interior 
points  to  take  place  before  any  contribution  was  registered. 
The  strains  at  the  integration  points  are  now  found  by 
taking  derivatives  of  the  assumed  displacement  functions. 


At  present  the  most  popular  method  for  Finite-Element 
elasto-plastic  analysis  is  the  incremental  Tangent  Modulus 
or  Variable  Stiffness  approach.  This  was  initiated  in  1965 
by  Pope  (46)  and  Swedlow  and  Yang  (47),  and  further  developed 
for  2-D  problems  by  Reyes  and  Deere  (48)  and  Marcal  and  King 
(49)  amongst  i.chers.  In  this  method  the  load  is  again  applied 
incrementally  but  now  it  is  the  matrix [D]in  equation  (1) 
which  is  updated  at  each  step  to  its  tangential  value,  [D-^]  » 
obtained  from  the  effective  stress-strain  curve  and  corresponding 
to  the  stress  state  reached.  This  means  that  a  new  set  of 
coefficients  for  the  equilibrium  equation  (ii)  must  be  found 
and  a  new  inversion  of  this  matrix  performed  at  each  load 
step.  A  diagram  of  the  process  is  given  as  Figure  lc). 

A  refinement  occasionally  employed  to  improve  accuracy  is 
described  by  Hcfmeister  et  al.  (50)  in  the  context  of  a 
2-D,  large  strain  analysis.  Here  an  equilibrium  check  is 
performed  at  each  load  step  (or  less  frequently)  and  thus  the 
nodal  point  equilibrium  error  reduced  by  using  a  Newton-Raphson 
interation  procedure.  A  comparison  between  Initial  Strain  and 
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Tangent  Modulus  methods  was  given  in  1968  by  Marcal  (51). 

He  demonstrated  that  a  larger  load  increment  could  be  taken 
for  a  variable  stiffness  approach,  but  as  the  modifications 
to  the  equilibrium  equations  required  at  each  step  are 
considerably  more  extensive  with  this  technique  the  overall 
efficiency  for  the  two  methods  remains  similar. 


Constant  stress  and  strain  elements  have  often  been 
used  .with  the  Tangent  Modulus  method  for  plane-stress  elasto- 
plastic  analyses.  Richard  and  Blacklock  (52)  reported  on 
such  an  element,  whereas  Felippa  (53)  described  work  using 
quadratic  functions  for  both  displacements.  A  linear  variation 
of  [Dr]  over  the  area  was  taken  with  values  calculated  at  the 
vertices  using  average  stresses  from  the  elements  at  each 
node.  The  problem  of  evaluating  the  matrix  [D-r]  for  the  Tangent 
Modulus  approach  is  equivalent  to  the  problem  of  computing 
the  pseudo— forces  for  the  Initial  Strain  approach.  Recently 
Bergen  and  Clough  (54)  considered  the  Felippa  approach  for  use 
with  a  quadrilateral  plate-bending  element,  but  concluded  it 
to  be  over-complicated  for  this  application.  In  addition, 
they  expressed  concern  about  assuming  a  continuous  |Dt1 
variation  when  this  is  frequently  not  so  within  the  element. 

The  approach  selected  by  Bergen  and  Clough  was  to  use  12 

subtriangles  per  element  and  to  take  values  of  stresses 

at  their  centroids.  Strip  integration  was  used  over  the  area 

and  Gaussian  quadrature  through  the  thickness.  As  with 

the  Initial  Strain  analysis  by  Levine  et  al. (45),  the  strains 

were  obtained  from  the  derivatives  of  the  assumed  displacement 

functions. 

A  good  deal  of  work  using  the  Tangent  Modulus  approach 
particularly  applied  to  shells  of  revolution  has  been 
performed  by  Marcal  (55,  56,  57,  58).  He  used  the  incremental 
theory  of  plasticity  and  the  von  Mises  yield  criterion. 
Isotropic  strain  hardening  was  assumed  and  a  linear  incremental 
elasto-plastic  stress-strain  relation  obtained.  Marcal  (58) 
evaluated  the  stresses  at  the  midway  point  of  his  axisymmetric 
shell  element  and  used  11  numerical  integration  points  through 
the  thickness.  In  1971  a  dynamic  analysis  of  beams  and  ring 
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was  presented  by  Wu  and  Witmer  (59).  They  assumed  an  elastic- 
perfectly  plastic  material  and  used  a  mechanical  sublayer 
model.  The  equations  of  motion  were  solved  by  a  time-wise, 
central-difference,  numerical  integration  ~>rocedure  using  a 
time-step  of  1  micro-second. 

A  further  Finite-Element  elasto-plastic  method  has  recently 
been  introduced  by  Zieniewicz  et  al.(60).  This  is  known  as 
the  incremental  Initial  Stress  approach,  in  which  plasticity 
is  accounted  for  by  adjustments  of  the  initial  stress  vector, 
tfg  ,  of  equation  (  4 .  The  approach  is  basically  a  modified 
Newton-Raphson  procedure  with  the  original  value  of  the 
stiffness  matrix  being  used  throughout.  Therefore  most  of 
the  advantages  of  the  Initial  Strain  approach  are  retained; 
indeed  the  two  methods  are  in  many  ways  parallel.  However, 
the  ability  to  analyse  perfectly  plastic  materials  is  a  notable 
additional  advantage  of  the  Initial  Stress  approach.  For 
both  these  procedures,  convergence  may  be  accelerated  by  over¬ 
adjustment  at  each  step.  Such  an  accelerator  for  use  with 
the  Initial  Stress  approach  has  been  developed  by  Nayak  and 
Zienkiewicz  (61).  Another  possibility  is  to  occasionally 
update  the  stiffness  matrix  to  its  tangential  value. 

Although  the  Initial  Stress  approach  has  been  demonstrated 
to  be  convergent  and  efficient  in  plasticity  problems  (60), 
it  may  not  be  as  suitable  for  problems  involving  large 
deformation  (62). 

Finally,  we  shall  mention  a  very  recent  elasto-plastic 
analysis  performed  by  Stricklin  et  al,(63).  in  which  Finite 
Element  and  Finite  Difference  methods  are  combined.  In 
this  formulation  the  effect  z£  the  non-linearities  are 
evaluated  by  means  of  Finite  Difference  expressions,  but 
the  general  set-up  of  the  analysis  is  Finite  Element.  The 
authors  present  results  for  shells  of  revolution  and  claim 
the  computational  procedure  to  be  an  order  of  magnitude 
faster  than  other  analyses. 


3.  MATERIAL  NON-LINEARITIES 


3.1  General  Theory  of  Plasticity 


Firstly  we  will  develop  a  general  theory  of  plasticity 
to  be  used  in  the  analysis.  A  yield  surface  can  be  defined 
by  the  general  equation 

f  (e,*)*o  (1) 

where  Tt  is  a  work  hardening  parameter.  Yielding  can  only 
occur  when  the  stress  components , J7  ,  satisfy  this  criterion 
for  the  instantaneous  value  of  %The  particular  yield  surface 
which  will  be  used  is  that  proposed  by  von  Mises.  For  a 
plate  or  shell  where  the  normal-to-plane  stress  component,  dz > 
is  negligible,  this  is  given  by 

+  <rf - ^ ^  4  3 +  ^4 )]*->( ft)  (2) 


where  Y(K)is  the  corresponding  value  of  the  uniaxial  yield 
stress. 


The  following  incremental  flow  rule  of  plasticity  is 
employed  in  the  analysis:- 


--  X 


(te) 


Here  ^represents  the  vector  of  incremental  plastic  strain 
components  and  X  is  an  unknown  proportionality  constant. 
Equation  (3)  is  often  referred  to  as  the  Normality  Principle 
as  it  may  be  interpreted  as  the  requirement  that  ^must  be 
normal  to  the  yield  surface  in  n-dimensional  stress  space. 


Explicit  formulations  of  the  incremental  plastic  stress- 
strain  relations  have  been  derived  by  Yamada  et  al.  (66)  and 
Zienkiewicz  et  al.  (60).  They  assume  that  for  an  infinitesimal 
stress  increment  the  corresponding  changes  in  strain  may 
be  divided  into  elastic  and  plastic  parts,  as 

r  ile.  Aj't*  t  a \ 


is  ■  m- is 
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By  differentiation  of  the  yield  surface  and  the  definition  of 

A  -  <AK  J- 

”  dll  X  (6> 


we  can  obtain  the  matrix  equation 


1 

2f 

1 

d<rK 

1 

9? 

! 

PE 

1 

dtx* ) 

2b 

I 

p'tz-x 

1 

2£ 

J 

^*3 

at 


dTxij 

X 


„  af  if'  af  if  4  \ 

[  0  J  [  d<r»  9ff$  n*K  7>tiy  n  J  [A  j 

fhis  type  of  relationship  has  been  used  with  success  by  Marcal 
and  King  (49)  in  their  Incremental  Tangent  Modulus  approach. 

For  our  purposes  however,  it  is  more  convenient  to  eliminate 
X  as  demonstrated  by  Zienkiewicz  et  al# (60).  to  give 


Jio~  —  [Plep  iS 


where 


[Dlr  [D] 


S  -  H 


-  X 


w(|(|)’[ol 


30T«(i 


(10) 


Here  H  is  the  slope  of  the  uniaxial  stress  versus  plastic 
strain  plot  at  the  value  of  the  initial  effective  stress. 

In  the  present  case,  using  the  von  Mises  yield  surface 
given  by  equation  (2)  we  have 


13 


. 


s 


where 

[Fm]  - 


If  we  define 


0  0  3 

0  0  0  s 

0  0  0  0  3 


^  =  [O]  [Fm]  £ 


equation  (10)  becomes 

5=  H'+/[F„]T2 


(11) 


(12) 


(13) 


(14) 


and  we  can  re-write  equation  (9)  as 

[Ply-  [0]  (is) 


3. 2  Finite  Element  Elasto-Plastic  Procedure 


The  method  chosen  to  account  for  plasticity  effects  is 
the  incremental  Initial  Stress  approach  suggested  by 
Zienkiewicz  et  al.  (60)  with  adaptions  to  give  compatibility 
with  both  the  hybrid  stiffness  formulation  and  dynamic  response 
analysis.  The  principal  advantage  of  this  formulation  is 
that  the  original  elastic  material  properties  are  used 
throughout  the  analysis  witn  plasticity  being  included  by 
means  of  a  system  of  pseudo-forces  acting  at  the  nodes. 


For  each  finite  element,  e,  the  subroutine  ELPL  is  called 
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by  the  master  programme  to  compute  these  forces,  JR  ,  from 
the  nodal  displacements  accumulated  during  a  time  increment. 
Gaussian  quadrature  < *  67)  is  used  to  obtain  Rl  from 
residual  (initial)  stress  values  computed  at  discrete  stations 
throughout  the  volume  of  the  element. 

In  the  original  Initial  Stress  procedure  proposed  for 
static  analysis,  iteration  we^s  performed  within  each  time 
increment  until  the  pseudo-forces  became  negligibly  small. 

This  is,  however,  believed  to  be  unnecessary  for  dynamic 
response  analysis  where  these  forces  may  be  carried  over 
to  act  during  the  next  time  increment. 

3.3.  The  Elasto-Plastic  Subroutine  ELPL 

A  flowchart  for  this  subroutine  is  included  in  Appendix 
I  and  the  operations  it  performs  will  now  be  described.  The 
process  is  represented  diagramatically  in  Fig.  2. 


Step  1 


From  the  nodal  displacements,  accumulated  for 

an  element  during  a  time  step  n ,  n  •*  I  the  change  in  the  stress 


coefficients  can  be  found  as 


1$  -  L< 


(16) 


Here  the  matrix  jcj  has  been  derived  during  the  normal  hybrid 
stiffness  computations. 


For  a  Gauss  point,  g,  use  the  hybrid  stress  assumptions 
to  obtain  elastic  increments  of  stress  and  strain, 


-  fc]3£f 


(17) 


h5  --  M  & 

^  (18) 

Add  dcr  to  the  existing  stress  components,  &  ,  to  give  O ?  . 

tm  15  — 


Step  3 

Test  whether 


p  (zl  ,*0  >  o 


(19) 


Herelt  refers  to  its  value  at  the  start  of  the  increment.  If 
expression  (19)  is  not  satisfied  then  the  increment  is 
entirely  elastic.  Therefore  up-date  stress  and  strain 
components  as 

JZ,  =  (70) 


(21) 


and  proceed  to  next  Gauss  station. 


If  expression  (19)  is  satisfied,  test  whether 


(22) 


using  the  stresses  at  the  start  of  the  increment.  If  expression 
(22)  is  also  satisfied  proceed  to  next  step;  if  not  interpolate 
for  proportions  of  elastic  stress  and  strain  occuring  above 
yield,  revise  Sv*  and  to  these  values  and  update  £? 

O 

and  6  to  the  yield  surface. 


Step  4 

t 

Call  the  Subroutine  HDASH  to  calculate  H  and  hence 
evaluate  the  scalor,  S,  from 


5=  H'4 


(23) 


and  the  elasto-plastic  matrix 


foj.  =  [o]  -  £ 


(24) 


Now  the  elasto-plastic  stress  increments  may  be  calculated 


from 


hV-  K  & 


(25) 
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As  equation  (25)  is  valid  for  infinitesimal  strain  increments 
only,  it  is  possible  that  for  a  finite  step  the  stresses 
calculated  wiil  slightly  exceed  the  yield  surface.  If 
this  occurs  the  stresses,  Sa^,}are  scaled  down  to  the  yield 
surface. 


Step  5 


Calculate  the  residual  (initial)  stress  vector 


~  <l4 


and  update  the  stress  components  as 


<7  -  £Te  -  L^a 


(27) 


and  the  strain  components  as  equation  (21), 


Step  6 


Repeat  Steps  2-5  for  each  Gauss  station,  hence  evaluating, 
by  Gaussian  quadrature,  the  nodal  pseudo-force  vector  in 
local  co-ordinates  as 


(28) 


The  derivation  of  this  expression  is  presented  in  Section  6. 


Step  7 


Transform  pseudo-forces  to  global  axes  by  performing 


r "  -  [li 


and  accumulating  to  form  ft 


4.  GEOMETRIC  NONLINEARITIES 


As  soon  as  the  maximum  transverse  displacement  of  the 
panel  reaches  about  one-half  the  panel  thickness  the  non¬ 
linearity  caused  by  the  large  displacements  should  be  accounted 


for.  An  incremental  process  originated  by  Argyris  et  al. 
(68,  69)  is  employed.  The  element  stiffness  matrices 
are  systematically  revised  throughout  the  analysis  by 
transforming  the  original  matrices  fjOc>  to  the  deformed 
orientations  using  the  equation 


(30 


This  transformation  accounts  for  gross  movement  of  the  whole 
element  as  a  rigid  body.  In  addition,  the  original  equilibrium 
of  each  element  is  disturbed  by  the  effects  of  straining 
and  large  rotations  from  the  initial  configuration,  and  to 
balance  this  a  vector  of  cumulative  nodal  pseudo-forces, 

S  ,  is  included  in  the  equations  of  motion  of  the  system. 


At  the  end  of  each  step,  the  forces  due  to  additional 

f  C  •'Ml 

straining,  0$e  may  be  evaluated  as 


fS 


n,H4i 


-  DC*  a* 


(31) 


Here 


L*T 


is  the  Kvised  stiffness  >.t  the  half-step  and  its 

significance  in  the  solution  procedure  is  explained  in  Section 

5.  In  addition,  forces  incorporating  large  rotation  effects 

and  acting  on  each  element  are  found  by  transforming  the 

r«>»% 

local  pseudo-force  vector,  , at  the  beginning  of  the  time 
step  as 


S&  «  S  ^ 


Ssl  =[t] 


Sl 


(32) 


Hence 


H4«  rfc-H  crK,n+’ 

S  =  S.  -v 


(33) 


Thus  it  is  now  possible  to  include  both  material  and  geometric 
non-linearities  into  the  analysis  without  entering  the 


necessarily  time-consuming  elemental  stiffness  routines  more 
than  once. 

5.  NUMERICAL  INTEGRATION  OF  THE  EQUATIONS  OF  MOTION 


The  equations  of  motion  of  the  system  can  be  written  as 


+[&3_|  +  [k]  ^  =  pw 


To  represent  the  non-linear  system  it  is  necessary  to 
adopt  a  step-by-step  procedure  taking  linear  conditions  to 
prevail  throughout  each  small  time  interval.  Thus  during 
a  particular  time-step  n,n+i  we  can  write. 


n,n*i 


(  35) 


The  scheme  selected  for  the  numerical  integration  is  an 
incremental  predictor-corrector  procedure  using  Runge-Kutta 
extrapolation  techniques  of  O(h^)  truncation  errror.  This 
formulation  conveniently  permits  the  inclusion  of  both  the 
material  and  geometric  non-linearities  present. 


For  transient  response  analysis  damping  has  very  little 
effect  on  the  general  solution  and  for  the  present  will  be 
excluded  from  the  analysis.  If,  however,  in  the  future  it 
is  desired  to  include  damping,  only  very  slight  modification 
of  the  programme  would  be  required. 

During  one  time-increment  the  following  scheme  of 
operations  is  performed: 


Step  1 

Compute  the  acceleration  components  vector  at  the 
beginning  of  the  time  step  as 


and  hence  evaluate  storage  vectors 

a*.  -  4"  h  +  v,” 


J, 


* 


(36) 


(37) 
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Step  2 

Extrapolate  for  incremental  displacements  at  the  half-time 
interval  as 

=•  i  U  ■+  -g  k  y*  (38) 


and  hence  compute  the  corresponding  accelerations  as 


Evaluate 


Va  =  ^ 


and  update  storage  vectors  as 

+  i  Ys 

and  *- 


(39) 

(40) 

(41) 

(42) 


Step  3 

Extrapolate  for  revised  incremental  displacements  at  the 
half-time  interval  as 


Hence  compute 
and 


Evaluate 


(43) 

C-  h? 

(44) 

(45) 

t  114*4 

v,3  =  k4* 

(46) 

and  update  storage  vector  to  its  final  incremental  value  as 

94* iyj)  (a?) 

and  update  as 

Aj,  =  ^  +  2  MJ 
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(48) 


Extrapolate  for  incremental  displacements  at  the  end 
of  the  time-step  as 


(49) 


hence  compute 


V*  MW  . 

h") 


(50) 


Evaluate 


..M4I 


|  •  •  T  I 

V4  =  k  4- 

^  T  (51) 

A  *’ 

and  update  storage  vector  /J^to  its  final  incremental  value  as 


£4*  *  tin* 


Step  5 


Update  the  displacement  vector 


retaining  in  store 

Update  velocity  vector 

*  #  h 


’  '  h  ,  Av  ' 

+£3;* 


(52) 


(53) 


(54) 


Step  6 

Compute  the  nodal  pseudo-forces  incorporating  geometric 
o  I 

non-linearities,  j  *  as  described  in  Section  4. 

Transform  the  incremental  nodal  displacements  to  local 
axes  for  each  element  using 


(55) 


and  hence  update  the  vector  of  pseudo-forces  for  each  element  as 

e,*  r  .  e. 


(56; 


Step  7 


For  each  element  call  the  subroutine  EL.eL  to  calculate 
the  nodal  pseudo-forces  accounting  for  plasticity,  , 
from  t  as  described  in  Section  3.3. 


6.  Derivation  of  Plasticity  Pseudo-Force  Expression 


It  has  been  demonstrated  by  Zienkiewicz  (.60) 
for  an  aelement,  the  vector  of  nodal  pseudo-forces 
residual  (initial)  stresses  can  be  expressed  as 

Rl4[B]T^3  Jv 

where  [bJ  is  the  matrix  relating  strain  and  nodal  displacements  as 


6 

1 — * 

(58) 

For  the  Hybrid 

approach  we  have 

4- 

■  H  V 

(59) 

cr 

(60) 

and 

e 

=  M  tQ]  [c]  X 

(61) 

that, 
due  to 

(57) 


Hence  by  comparison  with  equations  (58)  and  (57)  we  can  write 


(62) 
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7.  Mass  Matrix 


It  has  been  shown  by  Dungar  et  al.  (65)  that  in  order  to 
achieve  improved  results  in  dynamic  problems,  the  mass  and 
hybrid  stiffness  matrices  should  be  consistent  with  respect 
to  assumed  generalized  displacement  patterns.  In  the  hybrid 
stiffness  matrix  these  are  taken  along  the  element  boundaries 
only,  so  it  is  necessary,  when  formulating  the  mass  matrix, 
to  extend  the  functions  to  cover  the  whole  area  of  the  element, 


If  q.  and  u  are  respectively  the  vectors  of  generalized 
nodal  and  boundary  displacements  then  they  may  be  related  by 
a  matrix/ V  las  in  the  equation 

M  -Lv]^  (63) 

If,  also,  the  vector  of  generalized  displacement  within  the 
element,  d,  is  related  to  u  by: 


d  =  DO  a 


(64) 


then  the  element  m.ass  matrix  may  be  written  as: 


Dir  -  ?  [vj  d*  [v] 


(65) 


where 


p  =  mass  density  of  the  material  of  the  element. 
If  we  put  r 

[Ja]  =  f  [J]T[1]  d* 

then  A 


(66) 


in r=  f[v]T[Jc]fr] 


(67) 


The  degrees  of  freedom  at  each  node  follow  the  ordering 
w,  the  transverse  displacement  and  the  two  rotations'^:  and 
-$y  about  the  x  and  y  axes  respectively.  For  the  general 
triangular  element  orientated  with  respect  to  the  local  axes 
as  in  Fig.  3,  the  ,v]and  [jjp matrices  are  as  shown  in  Fig.  4 
and  5  respectively.  Transformation  to  global  axes  is  achieved 
by  the  same  process  as  the  stiffness  transformation. 


8.  BLAST-WAVE  REPRESENTATION 

The  subroutine  IMP  calculate-s  the  loading  vector,  r  » 

rv' 

at  time  increment,  n,  from  the  impulsive  forcing  function, 

P(t).  A  number  of  different  functions  have  been  tested  for 
their  ability  to  represent  the  pressure-time  curve  for  blast 
loading  (fig.  6).  Firstly,  polynomial  curve  fitting  was 
tried.  This  gave  low  root-mean-square  errors  for  the  higher 
order  polynomials,  but  the  representation  of  available  experimental 
data  gave  alternatively  over  and  under  estimates  along  the 
length  of  the  actual  curve.  Secondly,  an  exponential 
function  given  in  Ref.  70  was  tested: 

— 

?(t)  '  f>o  (  1  -  Vfcf )  e  (68) 

Here  is  the  initial  (maximum)  blast  pressure  and  t  is  the 
intercept  with  the  time  axis.  Fig.  3  shows  that  this  function 
overestimates  the  absolute  value  of  pressure  throughout. 

Finally,  a  modification  was  made  to  the  expression  (68)  to 

9ive  A  _  -fc/j l  A  -fc/. 


(i  -  Vtf)  -% 

P(fc) =  JV(|  +  t/tc)  € 


(69) 


This  function  is  seen  to  give  much  better  agreement  with 
the  actual  curve,  and  has  been  programmed  as  subroutine  IMP. 
When  values  for  p0  ,  t,  t/and  t*  are  input  the  subroutine  will 
calculate  the  value  of  the  loading  vector  at  the  time,  t. 

9.  RECOMMENDATION  FOR  FURTHER  WORK. 


Whilst  no  directly  applicable  results  have  been  obtained 
to  date,  the  authors  have  been  encouraged  by  the  rapid  develop¬ 
ment  of  the  various  subroutines,  see  Appendix  II.  The  algo- 
rithmns  employed  for  numerical  integration  are  well  known  for 
th^ir  numerical  stability,  and  the  progressive  nature  of  the 
spread  of  plasticity  through  each  element  should  improve  or 
other  elasto-plas tic  methods. 

It  is  recommended  that  the  work  is  continued  to  fully 
exploit  these  developments. 
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APPENDIX  I  -  FLOWCHARTS 


Call  TftfHVS  ^[L]0,hfn>K.a; 
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APENDIX  II  -  SUBROUTINES  TESTED  AS  AT  31.8. '72 


Subroutine  EQUIVSTRESS 

Gien  a  vector  of  stress  components,  at  a  Gauss  station 
this  subroutine  will  calculate  the  value  of  the  corresponding 
equivalent  stress,  (7?  as 


a  = 


(  A.  1 ) 


Subroutine  DEPS 


This  subroutine  calculates  the  incremental  value  of  the 
equivalent  plastic  strain at  a  Gauss  station,  based  on 
a  vector  of  corresponding  incremental  plastic  strain  components, 

Set. 


*  (t [  +  +  e  [  ( <5  v  , 


Subroutine  UNISTR 


(A.  2) 


The  equation  of  the  uniaxial  stress  (  Y(K)  )  versus 
uniaxial  plastic  strain  <  )  curve  is  taken  to  be 

«  A  (y(K) -*i.)g  (a. 3) 

where  is  the  uniaxial  stress  at  first  yield. 

When  given  values  for  A,  B7o*M,  and  subroutine  UNISTR 
will  calculate  the  value  of  the  uniaxial  stress,  .  This, 

together  with  the  equivalent  stress,  (J  ,  is  required  in  the 
evaluation  of  the  yield  surface  function  at  a  point,  as 


F(e,*0  =  c  -  YfJt) 


(A.  4) 


Subroutine  HD ASH 

i 

This  subroutine  calculates  the  value  of  H  ,  which  is 
the  instantaneous  slope  of  the  uniaxial  stress  versus  plastic 
strain  curve.  This  may  be  found  from  equation  (A. 3)  and  >5 
given  by  ( 

u'=  4 yfic) .  1  fep.»]s~ 

a  fl.E  \  /\  )  (A-5) 

Subroutine  MIN^ 


Two  matrix  inversion  routines  which  are  particularly 
suitable  for  the  present  problem  have  been  programmed  and 
tested.  The  more  accurate  of  the  two  empoys  full  pivotting, 
but  requires  more  storage  and  computation  time  that  the  other 
which  uses  partial  pivotting.  Both  routines  w j 1 X  be  tested 
in  the  overall  programme  and  the  better  one  selected  at  that 
time  when  the  accuracy  -  ’cost*  balance  can  be  assessed. 

Subroutine  TRANS 


Here  the  matrix  [L]  is  calculated  based  on  a  given  set 
of  nodal  co-ordinates.  This  matrix  performs  the  transformation 
from  global  to  local  co-ordinates  e.g. 


Vector  algebra  is  used  to  form  the  elements  of  [l])0. 


(A.  f>) 


Subroutine  KTRAN 


Using  the  matrix]^ computed  in  subroutine  TRANS,  this 
routine  sets  up  the  matrix  [Tj  given  by 

fT]0= 

p  L  ,  (A. 7) 


0  L 
0  0  L 
P  0  p  L 

for  a  four  noded  element . 
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KTRAN  then  performs  the  transformation  of  the  elemental 
stiffness  matrix [kJ  from  its  original  axes  to  global  axes,  as 


do:  =  ur  m:  [t] 


(A. 8) 


Subroutine  ELPL 

This  subroutine  computes  the  elemental,  plastic,  pseudo- 
force  vector,  R*.  ,  from  the  corresponding  incremental, 
nodal  displacement  vector  A  refined  flow  chart  for 

ELPL  is  presented  in  this  report. 


The  central  operation  in  this  subroutine  is  the  integration 
of  the  residual  (initial)  stresses,  Afloat  the  stations  to 
yield  Rt  .  For  each  station  the  amount  of  computation  and 
storage  is  great  and  hence  the  most  suitable  quadrature  scheme 
will  be  the  one  which  uses  the  minimum  number  of  stations  for 
a  given  accuracy.  The  obvious  choice,  therefore,  is  Gaussian 
quadrature,  and  when  we  apply  it  to  the  present  problem  we 
obtain: 

r:  -aT  />*  (M  PflcT)r  U 

where is  the  product  of  Gauss  weights  over  the  area  and 
thLoughout  the  thickness  at  station  and  .A  is  a  factor 
adjusting  the  limits  of  integration  as  those  appropriate 
to  the  particular  element  concerned. 

Originally  it  was  expected  the  triangular  elements  would 
be  the  most  suitable  for  use  in  the  analysis,  and  that  seven 
stations  over  the  area,  and  a  Gauss-Hammer  quadrature  scheme 
(71),  would  be  required.  However,  it  now  appears  that  by 
usipg  quadrilateral  elements  and  just  four  stations  over  the 
are.'i  comparable  results  can  be  achieved  with  much 
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less  storage  and  computation  time.  This  is  primarily  because 
the  stiffness  matrix  ,  [K]  ,  for  a  four-noded  element  is 

appreciably  more  accurate  than  that  for  a  triangle,  and  it 
therefore  follows  that  the  same  is  true  for  the  residual 
stress  vector,  Ac,  at  the  chosen  number  of  Gauss  stations. 

Stress  variations  through  the  thickness  of  the  elasto- 
plastic  element  are  generally  more  complex  than  those  over 
the  area,  and  hence  in  the  first  tests  seven  stations  through 
the  thickness  are  being  employed.  As  the  first  yielding  of 
a  panel  subjected  to  lateral  loading  occurs  on  the  surface, 
it  is  desirable  that  the  extreme  stations  through  tie  thickness 
are  near  to  the  surface  in  order  that  the  initial  yielding 
is  quickly  detected.  By  using  seven-point  Gaussian  quadrature 
the  distance  of  the  outer-most  stations  from  the  surface  is 
just  2.5%  of  the  plate  thickness. 

Subroutine  KMAT 


This  subroutine  calculates  the  elastic  stiffness  of  a 
triangular  or  quadrilateral  shell  element  with  respect  to 
a  global  co-ordinate  system.  The  method  used  is  known  as 
the  "Hybrid"  technique,  and  the  element  has  been  well  tested 
in  plate  and  shell  situations. 

Subroutine  MMAT 


This  calculates  the  mass  matrix  for  a  general 
triangular  element  and  transforms  it  to  global  axes.  The 
theory  is  given  in  Section  7. 
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Fig  J-  Diaqrci  nr\  m  a  •  i  c  Repress  y\fat\£>r\  of  Imfial 
Stress'  E  taste- Phytic  Process  for  Jncre.ry.ekyt  h,  h-*  1 
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Fig.  4,  fv] for  General  Triangular  Element. 
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